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1.  Introduction 

As  a  result  of  their  high  hardness,  high  stiffness,  and  rela¬ 
tively  low  mass  density,  ceramic  materials  are  used  frequently  in 
industrial  and  military  applications.  Examples  of  the  latter  include 
ballistic  protection  systems,  i.e.,  armor,  for  vehicular  and  personal 
protection. 

A  ceramic  material  of  high  current  interest  is  boron  carbide,  B4C. 
Single  crystals  of  boron  carbide  exhibit  a  rhombohedral  (i.e.,  trig¬ 
onal)  crystal  structure:  space  group  R3m,  centrosymmetric  point 
group  3m,  and  Laue  group  Rl.  Various  polytypes  (i.e.,  structural 
variants)  exist.  In  the  boron  carbide  unit  cell,  C-C-C  or  C-B-C 
atomic  chains  are  aligned  with  the  [00  01]  direction  (hexagonal 
notation)  and  thread  the  basal  layers  of  icosahedra  located  at  the 
rhombohedral  vertices  (Yan  et  al„  2009).  Strictly,  the  rhombohedral 
space  group  description  applies  for  the  particular  polytype  with 
ideal  structure  B12(CCC),  denoted  herein  simply  as  “CCC”.  Slight 
distortion  from  a  trigonal  structure  may  occur  in  other  polytypes. 
The  polar  CBC  polytype,  with  structure  formula  B-nCp(CBC)  and 
denoted  herein  simply  as  “polar”,  is  most  thermodynamically  sta¬ 
ble  (i.e.,  lowest  ground  state  energy)  and  thought  most  abundant 
(Fanchini  et  al.,  2006);  stoichiometry  of  boron  carbide  may  also 
depart  from  the  ideal  boron-to-carbon  ratio  of  B4C,  with  lower 
carbon  content  possible  (Dodd  et  al.,  2002). 
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Polycrystalline  boron  carbide,  with  grains  of  typical  sizes  on 
the  order  of  10p,m,  has  a  low  density  of  2.52  g/cm3,  high  elastic 
modulus  of  470  GPa  (McClellan  et  al.,  2001),  high  hardness  on  the 
order  of  30  GPa,  and  a  large  Hugoniot  elastic  limit  (HEL)  on  the 
order  of  15-20  GPa  (Vogler  et  al.,  2004;  Zhang  et  al.,  2006).  Duc¬ 
tility  of  boron  carbide  is  low;  though  partial  dislocations,  stacking 
faults,  and  twins  have  been  observed,  such  defects  are  thought  to 
be  relatively  scarce  and/or  immobile  during  high  rate  deformation 
in  commercial  material  produced  by  hot  pressing.  Instead,  under 
impact  loading,  strength  and  stiffness  of  the  material  are  thought 
to  be  limited  by  amorphization  and  cleavage  fracture  (Chen  et  al., 
2003).  Shock-induced  amorphization  of  boron  carbide  was  dis¬ 
covered  during  post-mortem  examination  of  fragments  recovered 
from  ballistic  experiments  (Chen  et  al.,  2003),  where  it  was  sug¬ 
gested  that  onset  of  amorphization  may  be  associated  with  loss  of 
intrinsic  stability  in  the  sense  of  Born  (1940).  Previously,  association 
of  the  Born  instability  with  amorphization  was  posited  for  quartz 
(Gregoryanz  et  al.,  2000),  which  like  boron  carbide,  has  a  rhom¬ 
bohedral  crystal  structure  and  demonstrates  a  decrease  in  elastic 
shear  stiffness  C44  at  high  pressure. 

Density  Functional  Theory  (DFT)  calculations  of  hydrostatic 
loading  (Fanchini  et  al„  2006)  predicted  collapse  of  the  crystal 
structure  of  boron  carbide  into  amorphous  regions  consisting  of 
a  mixture  of  carbon  and  boron  icosahedra,  with  segregation  most 
energetically  favorable  for  the  CCC  polytype.  First  principles  cal¬ 
culations  (Yan  et  al.,  2009)  implied  that  amorphization  follows 
cross-linking  of  the  C-B-C  atomic  chain  of  the  boron  carbide  unit 
cell  with  icosahedral  atom(s)  at  a  vertex  of  the  rhombohedral  unit 
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cell  in  the  polar  polytype.  The  cross-linked  structure  is  unstable  at 
ambient  pressure;  hence,  amorphization  and  fracture  may  occur  on 
decompression.  Cross-linking  occurred  for  uniaxial  straining  along 
the  c-axis  ([0  00  l]-direction)  at  an  axial  pressure  around  20GPa 
but  was  not  observed  for  hydrostatic  compression  to  60  GPa.  Amor¬ 
phization  has  also  been  observed  in  dynamic  indentation  (Subhash 
et  al„  2008).  Recent  DFT  results  (Taylor  et  al„  2011)  report  varia¬ 
tions  of  elastic  stiffness  components  of  several  polytypes  of  boron 
carbide  with  increasing  hydrostatic  pressure.  A  decreasing  tan¬ 
gent  elastic  coefficient  C44  with  increasing  pressure  was  predicted, 
leading  to  a  vanishing  eigenvalue  of  6  x  6  stiffness  matrix  C  at 
hydrostatic  pressure  p  ss  44  GPa  in  the  CCC  polytype  and  p  =»  68  GPa 
in  the  polar  polytype.  These  vanishing  eigenvalues  correspond  to 
onset  of  the  aforementioned  Born  instability. 

In  previous  work  (Clayton,  201 2),  a  third-order  anisotropic,  non¬ 
linear  elastic  continuum  model  was  developed  for  CCC  and  polar 
CBC  polytypes  of  boron  carbide  single  crystals,  using  a  combina¬ 
tion  of  published  data  from  experiments  and  first  principles  (Lee 
et  al„  1992;  McClellan  et  al.,  2001;  Yan  et  al„  2009;  Taylor  et  al„ 

2011) .  This  model  was  used  to  predict  the  onset  of  mechanical 
instability,  according  to  several  intrinsic  criteria  (Hill,  1975;  Parry, 
1978, 1980a;  Milstein  and  Hill,  1979;  Wang  et  al.,  1995;  Morris  and 
Krenn,  2000),  in  single  and  polycrystals  deformed  homogeneously 
and  elasto-statically.  It  was  determined  that  the  particular  intrinsic 
“B-criterion”  corresponding  to  incremental  loading  by  the  Cauchy 
stress  predicted  the  earliest  onset  of  instability  and  correlated  with 
the  experimental  HEL  for  uniaxial-strained  polycrystals  (Clayton, 

2012) .  As  discussed  in  Appendix  A,  in  contrast  to  classical  elastic 
stability  criteria  that  depend  on  the  load  environment  and  may 
require  consideration  of  global  inequalities  (Ogden,  1984),  “intrin¬ 
sic”  or  “internal”  stability  criteria  (Parry,  1978,  1980a;  Morris  and 
Krenn,  2000;  Clayton  and  Bliss,  unpublished)  depend  only  on  the 
local  state  (tangent  stiffness  and  local  stress).  For  homogeneous 
stress  fields,  intrinsic  stability  may  be  sufficient  for  classical  sta¬ 
bility  if  conjugate  stress-deformation  measures  of  the  classical  and 
intrinsic  criteria  coincide. 

In  this  paper,  the  elasticity  model  of  Clayton  (2012)  is  extended 
to  explicitly  account  for  the  amorphous  phase.  In  Section  2,  a 
general  finite  strain  framework  is  first  developed  in  the  context 
of  continuum  thermomechanics  and  is  then  specialized  to  boron 
carbide.  Amorphization  is  modeled  to  occur  upon  attainment  of 
intrinsic  instability.  Discussed  in  Section  3  are  model  implemen¬ 
tation  in  a  dynamic  finite  element  (FE)  context  and  “mesoscale” 
simulations  of  shock  and  dynamic  uniaxial  compression  of  B4C 
polycrystals.  Conclusions  follow  in  Section  4. 

2.  Theory 


parameter: 

f(X,  t)  =  0VX  e  crystal,  f(X,  t)  =  1VX  e  amorphous; 

Ff  =  F?[f(X,  t)[,  Ff(0)  =  l;  (2.2) 

with  1  the  unit  tensor.  Helmoltz  free  energy  per  unit  reference 
volume  is  of  the  functional  form 

<F  =  «F[E(F£),  (?,£].  (2.3) 

Temperature  is  9\  Green  elastic  strain  is  E  =  (l/2)(  f£T  F£-l). 
Dependence  on  X  is  implicit  in  (2.3)  for  a  heterogeneous  polycrys¬ 
tal  whose  elastic  moduli  vary  from  grain  to  grain.  Interfacial  energy 
associated  with  phase  boundaries  is  omitted  in  this  description  that 
treats  £  as  an  internal  state  variable;  such  surface  energy  could  be 
included  by  extending  (2.3)  to  include  dependence  on  Vf,  as  in 
phase  field  theory  (Clayton  and  Knap,  2011a). 

Let  o',  P,  and  S  denote  Cauchy,  first  Piola-Kirchhoff,  and  second- 
Piola-Kirchhoff  stresses,  related  by  a  =J_1  PFT  =J_1  F5FT  with/  =  det 
F.  Let  l)  =  /  denote  particle  velocity,  Vt >  =  F,  and  U=>P  +  9r]  be 
internal  energy  with  r/  entropy  per  unit  reference  volume.  Let  po 
and  p  denote  reference  and  spatial  mass  densities,  and  Q  the  ref¬ 
erential  heat  flux  vector.  The  usual  balance  laws  of  continuum  field 
theory  (Clayton,  2011a)  apply,  expressed  here  in  the  absence  of 
body  forces  and  heat  sources  other  than  conduction: 

p0  =  pj,  W  P  =  p0V,  FPt  =  PFt,  L/  =  P  :  F  -  V  ■  Q.  (2.4) 


The  Clausius-Duhem  inequality  can  be  written 
P  :  F  -  —  9rj  -  0_1  Q  ■  V0  >  0.  (2.5) 


Use  of  (2.1)-(2.4)  in  (2.5)  and  applying  standard  thermodynamic 
arguments  (Clayton,  2011a)  results  in  the  hyperelastic  stress-strain 
law,  the  entropy  relation,  and  a  dissipation  inequality: 


9^  f-T 

9f£  9e 


9<F 


(F£TP)  :  ^  ^ 


9? 


f__Q.V0>O. 


(2.6) 


The  quantity  in  square  braces  can  be  viewed  as  a  thermodynamic 
driving  force  for  transformation. 

Let  a  superscript  0  enclosed  in  parentheses  signify  a  quantity 
for  a  pure  crystalline  phase  (£  =  0)  and  let  a  superscript  1  signify  a 
quantity  for  the  completely  amorphized  phase  (f  =  1 ).  Free  energies 
of  crystal  and  amorphous  phases  can  be  expanded  as  Taylor  series 
about  a  reference  state  wherein  (  £,  0)  =  (O,  90)\ 


A  general  continuum  constitutive  theory  is  developed  for  solids 
undergoing  a  crystal-to-amorphous  transformation  in  Section  2.1. 
This  theory  is  specialized  to  boron  carbide  in  Section  2.2. 

2.3.  Kinematics  and  thermodynamics 


«*0)  =  ne,  9, 0)  =  <p'0)  +  2,  c$LE„EKL 

+  JTCiwlmneUeklEmn  +  ■  ■  ■  +  S(0)(E,  9),  (2.7) 


Let  x  denote  the  spatial  position  of  a  material  particle  with  ref¬ 
erence  location  X.  Let  x  =  /(  X,  t)  denote  the  one-to-one  motion, 
with  t  time.  Letting  V  denote  the  material  gradient,  deformation 
gradient  F  is 

F  =  V/  =  F£Ff,  (2.1) 

where  £f  accounts  for  deformation  due  to  transforma- 
tion/amorphization  and  F£  is  the  remaining  elastic  deformation 
associated  with  local  stress  and  any  rigid  body  rotation.  Let  £  e  [0, 
1]  denote  an  order  parameter  distinguishing  between  crystalline 
and  amorphous  phases,  and  let  F^  depend  locally  only  on  this 


<F(1)=/-1«F(£,  9, 1)  =  V™ 

+  ^mEUE«L  +  3 ,  ^LMNEUEKLEMN  +  ■■■  +  S(1)(£,  9).  (2.8) 

Ground  state  energy  of  each  phase  (elastically  unstrained  and 
measured  at  9q)  is  <Fo-  Isothermal  second-  and  third-order  elas¬ 
tic  constants  (measured  at  9q)  are  Cyia  and  Cijklmn ■  Functions  S 
account  for  specific  heat  and  thermal  expansion  in  each  phase  and 
vanish  at  9  =  9g.  Energy  is  divided  by  J^  =  det  £f  in  (2.8)  so  that 
material  coefficients  are  defined  per  unit  unstressed  volume  of  the 
amorphous  phase. 
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Table  1 

Independent  nonzero  second-  and  third-order  elastic  constants  for  B4C  single  crys¬ 
tals  (GPa). 


Constant 

Polar  polytype 

CCC  polytype 

Reference 

C„ 

554 

487 

Taylor  et  al.  (2011) 

Cl2 

121 

117 

Taylor  et  al.  (2011) 

Cl3 

65 

66 

Taylor  et  al.  (2011) 

C33 

526 

525 

Taylor  et  al.  (2011) 

C44 

155 

133 

Taylor  et  al.  (2011) 

C14 

-18 

-18 

Lee  et  al.  (1992) 

Ko 

236 

222 

Clayton  (2012) 

Cm 

-873 

-667 

Clayton  (2012) 

Cl12 

-1668 

-1400 

Clayton  (2012) 

Cl13 

-974 

-753 

Clayton  (2012) 

Cl  23 

-3605 

-3913 

Clayton  (2012) 

Cl33 

-873 

-667 

Clayton  (2012) 

C222 

-873 

-667 

Clayton  (2012) 

C333 

-873 

-667 

Clayton  (2012) 

C344 

1645 

3015 

Clayton  (2012) 

2.2.  Model  and  parameters  for  boron  carbide  phases 

For  the  crystalline  phase  (£  =  0),  the  third-order  anisotropic  elas¬ 
tic  model  of  Clayton  (2012)  is  implemented.  Polar  (CBC)  and  CCC 
polytypes  are  considered  with  rhombohedral  symmetry.  Second- 
and  third-order  elastic  constants  in  the  undeformed  state  -  Cap 
and  Capy  in  Voigt  notation  -  are  listed  in  Table  1,  along  with  bulk 
modulus  1<0.  In  calculations  that  follow  in  Section  3,  temperature 
effects  on  the  mechanical  response  of  the  material  are  ignored,  cor¬ 
responding  to  3  =  constant  in  (2.7)  and  (2.8);  0  is  neither  explicitly 
prescribed  nor  calculated. 

Stress-induced  amorphization  in  B4C  is  not  well  understood,  and 
associated  physical  properties  are  not  fully  known.  A  number  of 
amorphous  structures  have  been  reported  from  DFT  (Fanchini  et  al„ 
2006;  Ivanshchenko  et  al„  2009;  Yan  et  al.,  2009)  with  disordered 
arrangements  featuring  various  atomic  clusters,  layers,  and  differ¬ 
ent  coordination  among  C  atoms.  Following  Ivanshchenko  et  al. 
(2009),  here  it  is  assumed  that  amorphous  and  crystalline  phases 
have  the  same  mass  density  po  =  2.52  g/cm3,  and  that  the  bulk  mod¬ 
ulus  of  the  glassy  phase  1<t  =  0.96I<0,  as  listed  in  Table  2.  Regarding 
the  former  assumption  (which  implies  ft  =  1 ),  other  DFT  work  (Yan 
et  ah,  2009)  suggests  that  mass  density  may  increase  on  structure 
collapse,  but  this  possibility  is  omitted  here.  In  the  absence  of  fur¬ 
ther  data,  it  is  assumed  that  f  £  =  1,  and  that  the  amorphous  phase 
behaves  as  an  isotropic  St.  Venant-Kirchhoff  solid  obeying  Cauchy’s 
symmetry  relations  (Clayton,  2011a)  leading  to  Poisson  ratio  v  =  1  /4 
and  shear  modulus  /r  =  (3/5)K!.  The  isotropic  assumption  is  valid 
for  amorphous  solids;  Cauchy  symmetry,  while  unproven  for  this 
material,  corresponds  to  pairwise  central  force  interactions  among 
atoms  and  is  characteristic  of  ideal  amorphous/glassy  materials 
(Pineda,  2006). 

With  these  model  assumptions  applied  to  boron  carbide,  free 
energies  (2.7)  and  (2.8)  become 

«*0)  =  ^o0)  +  l^pEaEp  +  ±Ca/,yEaEpEy, 

^1)  =  i^1)+^[(tr£)2  +  2£:£].  (2.9) 

Since  time  scale(s)  involved  in  amorphization  kinetics  are 
unknown  and  unavailable  from  static  DFT,  a  discrete  represen¬ 
tation  is  used,  and  the  domain  0  <  £  <  1  is  ignored  in  numerical 

Table  2 

Isotropic  elastic  properties  for  amorphous  phase. 

I<i  /K0  v  Reference 

0.96  \  Ivanshchenko  et  al.  (2009) 


simulations  reported  later.  Recrystallization  has  not  been  reported 
in  the  literature  and  is  precluded  by  the  model,  meaning  f :  0  ->•  1  is 
irreversible.  Following  earlier  conjecture  (Taylor  et  al.,  2011)  and 
analysis  (Clayton,  2012),  stress-induced  amorphization  in  B4C  at  a 
material  point  X  at  time  t  is  prescribed  to  occur  upon  attainment 
of  the  B-instability  criterion; 

det[B(X,  £)]  <  0=>  f(X,  t)->  1,  (2.10) 

where  incremental  modulus  B  is  defined  in  (A.7)  of  Appendix  A 
and  can  be  written  as  a  symmetric  6x6  matrix  (Wang  et  al.,  1995; 
Morris  and  Krenn,  2000;  Clayton,  2012).  Since  B  depends  on  the 
local  stress-deformation  state,  amorphization  propagates  at  the 
elastic  wave  speed.  This  transformation  model  is  advantageous 
because  no  fitting  parameters  are  introduced.  However,  in  this 
discrete  representation,  f  oo  at  the  instant  of  amorphization, 
precluding  pragmatic  enforcement  of  the  dissipation  inequality  in 
(2.6)  that  requires  time  differentiability  of  state  variables. 

3.  Numerical  method  and  results 

3.1.  Finite  element  method 

The  constitutive  model  of  Section  2.2  is  implemented  in  a 
dynamic  Lagrangian  finite  element  program.  Behavior  of  a  repre¬ 
sentative  microstructure  consisting  of  50  polyhedral  grains  (Fig.  1) 
is  studied.  The  FE  mesh,  with  sufficient  refinement  near  grain 
boundaries  and  triple  points  to  resolve  local  fields,  consists  of  1.6M 
hexahedral  elements.  Software  and  mesh  generation  procedures 
have  been  discussed  in  detail  elsewhere  (Clayton  et  al.,  2012a, b) 
in  the  context  of  other  ceramic  materials.  Faces  of  the  polycrys¬ 
talline  cube  are  normal  to  global  X,  Y,  and  Z  directions.  Four  initial 
microstructures  are  considered  with  random  lattice  orientation 
arrangements  labeled  1 , 2, 3,  and  4,  and  with  (  =  0  V  X  at  the  begin¬ 
ning  of  each  simulation.  Lattices  1  and  2  contain  mixtures  of  grains 
of  CCC  and  polar  poly  types ;  lattices  3  and  4  consist  only  of  randomly 
oriented  crystals  of  the  polar  polytype.  In  dynamic  simulations, 
transformation  criterion  (2.10)  is  enacted  at  an  integration  point 
when  the  minimum  eigenvalue  of  matrix  [B]  first  becomes  non¬ 
positive. 

3.2.  Shock  compression 

Boundary  conditions  establish  a  planar  shock  wave  that  tran¬ 
sits  across  the  sample.  The  edge  length  of  the  polycrystalline  cube 
is  40  p,m,  yielding  an  average  grain  size  consistent  with  the  real 
material  (Dandekar,  2001 ).  One  face  of  the  cube  is  assigned  a  uni¬ 
form  pressure  boundary  condition,  ramped  from  zero  to  p  over  a 
small  (4ps)  time  interval  and  then  held  constant  for  the  remain¬ 
der  of  the  simulation.  Lateral  boundaries  are  fixed  in  the  normal 
direction  to  maintain  a  uniaxial  strain  state,  while  the  face  oppo¬ 
site  the  pressurized  face  remains  free.  A  number  of  simulations  are 
conducted  with  p  of  10,  15,  20,  30,  40,  and  50  GPa,  applied  to  any 


(a)  (b) 


Fig.  1.  Synthetic  boron  carbide  polycrystal:  (a)  grain  structure,  transparent  view 
and  (b)  FE  mesh. 
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ofX,  Y,  or  Z  faces,  with  any  of  lattice  arrangements  1-4.  Note  that 
the  overall  stress  state  is  not  hydrostatic  since  p,  the  compressive 
normal  component  of  Cauchy  stress,  is  applied  to  only  one  face  of 
the  cube. 

Representative  results  for  a  polycrystal  consisting  of  both  poly¬ 
types  are  shown  in  Fig.  2;  analogous  results  for  a  polycrystal 
consisting  only  of  the  polar  polytype  are  shown  in  Fig.  3.  In  each 
case,  the  shock  wave  is  propagating  from  top  to  bottom,  and  has 
yet  to  transit  the  entire  sample  at  2  ns.  The  axial  component  of 
Cauchy  stress,  positive  in  compression,  is  shown  in  part  (a)  of  each 
figure,  and  can  exceed  the  applied  axial  pressure  p  by  f»15GPa. 
The  square  root  of  the  second  invariant  of  the  deviatoric  Cauchy 
stress,  s/h,  a  measure  of  the  magnitude  of  local  shear  stress,  has 
a  typical  value  on  the  order  of  lOGPa  and  can  locally  surpass 
15GPa  or  exceed  the  theoretical  strength  (sa(  1/1 0)C44 ).  Transfor¬ 
mation  in  the  mixture  of  CCC  + polar  polytypes  in  Fig.  2(c)  tends 
to  involve  entire  grains  (mainly  CCC),  whereas  transformation  in 
the  polar  polytype  in  Fig.  3(c)  is  less  pronounced  and  correlates 
with  areas  of  concentrated  axial  stress  (Fig.  3(a))  at  several  grains 
and  grain  boundaries.  Axial  stress  versus  volume  for  all  simula¬ 
tions  are  compared  with  experimental  data  (Vogler  et  al.,  2004; 
Zhang  et  al.,  2006)  in  Fig.  4(a).  Hugoniot  equations  (Clayton,  2010; 


Thurston,  1974)  have  been  used  to  estimate  “average”  volume  V  in 
the  shocked  state  and  the  “average”  shock  velocity  D  from  simula¬ 
tion  data: 


V  = 


(3.1) 


with  V0  =  (40  p,m)3.  Effective  particle  velocity  v  and  shock  stress  P 
are  estimated  from  simulation  data  as 

P^p,  (3.2) 

with  (Uf)  the  arithmetic  mean  of  particle  velocities  of  nodes  on 
the  free  face  of  the  polycrystalline  cube  at  t  =  5ns.  Recall  that  for 
uniaxial  strain  in  theXi  direction,  axial  Cauchy  stress  (i.e.,  the  Hugo¬ 
niot  “pressure”)  and  axial  first  Piola-Kirchhoff  stress  are  equal: 
On  =J~^Fu<Pv<=J-1FiiPu  =Pn.  where  Pjj  is  equivalent  in  magni¬ 
tude  to  P  in  the  Hugoniot  relations  (Thurston,  1974).  The  standard 
Hugoniot  relations,  with  the  material  in  front  of  the  shock  at  rest, 
are  derived  for  a  homogeneous  medium  with  a  planar  shock  (step 
discontinuity  in  stress)  in  one  spatial  dimension;  Eqs.  (3.1 )  and  (3.2) 
are  approximations  in  the  present  context  because  of  local  hetero¬ 
geneity,  the  small  but  finite  thickness  of  the  shock,  and  stress  wave 


Fig.  3.  Shock  compression  of  polar  polytype  (p  =  20GPa)  at  t=2  ns:  (a)  axial  stress,  (b)  shear  stress  and  (c)  order  parameter. 


vi  v0  v 

(a)  (b) 


Fig.  4.  Axial  or  shock  stress  ( p  or  P)  versus  axial  compression  ratio  V/Vo :  (a)  shock  simulations  (see  Section  3.2),  plate  impact  experiments  (Vogler  et  al.,  2004;  Zhang  et  al., 
2006)  and  (b)  uniform  compression  simulations  (see  Section  3.3). 
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Table  3 

Results  summary  at  r=5ns,  various  loading  directions  and  initial  microstructures, 
P  ~  p  =  20  GPa  shock. 


Load  dir. 

Lattice 

Initial  polar  frac. 

f 

v  (m/s) 

D  (km/s) 

VI  Vo 

X 

1 

0.45 

0.39 

604 

13.1 

0.954 

Y 

1 

0.45 

0.38 

622 

12.8 

0.951 

Y 

2 

0.77 

0.24 

626 

12.7 

0.951 

Z 

1 

0.45 

0.45 

635 

12.5 

0.949 

X 

3 

1.00 

0.09 

596 

13.3 

0.955 

Y 

3 

1.00 

0.03 

607 

13.1 

0.954 

Y 

4 

1.00 

0.05 

586 

13.5 

0.957 

Z 

3 

1.00 

0.06 

608 

13.1 

0.953 

interactions  with  boundaries  of  the  domain.  By  t=5ns,  the  main 
compressive  (shock)  wave  has  fully  interacted  with  the  free  surface 
(necessitating  the  factor  of  1  /2  in  (3.2)),  and  a  nearly  steady  value  of 
V  has  been  reached  in  each  simulation.  At  that  time,  some  unload¬ 
ing  of  the  sample  has  occurred  upon  reflection  of  the  compressive 
wave  at  the  free  surface,  but  results  reveal  that  further  transfor¬ 
mation  is  minimal  in  unloaded  regions  of  material.  This  behavior 
is  expected  since  transformation  in  B4C  is  driven  by  compressive 
stress  (Clayton,  2012). 

Experimental  values  for  Hugoniot  stress  P  and  volume  ratio 
VIVQ  =  pol P  are  obtained  directly  from  Vogler  et  al.  (2004)  and 
Zhang  et  al.  (2006).  In  these  experiments,  a  two-wave  structure 
was  observed,  whereas  in  the  present  simulations,  a  single-wave 
structure  was  observed.  In  simulations,  the  length  of  the  domain  in 
the  direction  of  shock  propagation  was  insufficient  to  permit  a  fully 
formed  sharp  shock  with  a  steady  profile  prior  to  arrival  at  the  free 
surface.  It  has  been  shown  elsewhere  (Clifton,  1971)  that  even  for 
an  isotropic  homogeneous  nonlinear  elastic-plastic  material,  the 
propagation  distance  needed  for  a  steady  wave  profile  to  emerge 
for  weak  shocks  can  exceed  0.5-1  mm.  Though  the  first  arrival  time 
of  the  elastic  wave  could  be  obtained  from  simulations,  it  was  not 
possible  to  extract  a  shock  velocity  D  (meaningful  for  comparison 
among  different  cases)  directly  due  to  the  small  but  finite  rise  time 
of  the  free  surface  velocity  and  heterogeneities  induced  by  grain 
structure  and  anisotropy. 

Comparison  with  experiment  is  reasonable;  at  each  impact  pres¬ 
sure,  V/Vo  for  the  polar  polytype  alone  exceeds  that  for  the  mixture, 
demonstrating  higher  average  stiffness  of  the  former  (i.e.,  amor- 
phization  tends  to  decrease  axial  stiffness).  Quantities  entering 
(3.1)  are  compared  for  all  20  GPa  simulations  in  Table  3;  £  is  the  vol¬ 
ume  fraction  of  amorphized  material  in  each  polycrystal  at  t=  5  ns. 
Inclusion  of  the  less  intrinsically  stable  CCC  polytype  correlates 
with  increased  transformation  and  reduced  shock  velocity.  The 
solid  curve  in  Fig.  4(a)  corresponds  to  a  uniform  compression  sim¬ 
ulation  as  discussed  in  Section  3.3. 

3.3.  Uniform  compression 

Boundary  and  initial  conditions  establish  an  average  uniaxial 
strain  state  with  a  nominal  strain  rate  of  105/s  (Clayton  et  al„ 
201 2a, b).  One  face  of  the  cube  is  assigned  a  constant  velocity  bound¬ 
ary  condition  inducing  axial  compressive  strain.  Lateral  boundaries 
are  constrained  in  transverse  directions  to  maintain  a  uniaxial 
strain  state,  while  the  face  opposite  the  moving  face  is  rigidly  fixed. 
An  initial  velocity  gradient  is  prescribed  consistent  with  the  applied 
strain  rate,  so  that  the  material  undergoes  shocldess  compression. 
A  number  of  simulations  are  conducted  with  strain  applied  in  any 
of  X,  V,  or  Z  directions,  with  any  of  lattice  arrangements  1-4.  The 
Courant  condition  restricts  the  integration  timestep  At  (propor¬ 
tional  to  the  minimum  element  dimension)  and  thus  the  maximum 
average  strain  that  can  be  attained  at  a  fixed  strain  rate  e.  For  exam¬ 
ple,  the  number  of  integration  cycles  (proportional  to  wall-clock 
compute  time)  to  achieve  a  strain  e  is  n  =  e/(eA t).  Therefore,  for 


uniform  compression  the  edge  length  of  the  polycrystal  cube  is  set 
to  1  mm. 

Characteristic  results  for  polycrystals  of  CCC  +  polar  and  polar 
polytype  only  are  shown  respectively  in  Figs.  5  and  6.  Local  con¬ 
centrations  of  axial  and  shear  stress  are  evident  in  certain  crystals 
and  grain/phase  boundaries.  Transformation  in  the  polar  polytype 
is  visibly  restricted  to  a  few  small  regions  (Fig.  6(c))  that  correlate 
with  stress  concentrations  at  grain  boundaries  in  Fig.  6(a)  and  (b). 
As  was  the  case  for  shock  loading,  transformation  in  the  CCC  +  polar 
mixture  (Fig.  5(c))  involves  comparatively  large  volumes  often  cor¬ 
responding  to  whole  grains.  Average  axial  stress  p  (normal  force  per 
unit  equivalent  reference  or  current  area)  acting  on  the  moving  face 
of  the  polycrystal  cube  is  compared  for  six  uniform  compression 
simulations  in  Fig.  4(b).  Undulations  correlate  with  amorphization 
events.  Despite  differences  in  absolute  size  of  specimens,  results 
in  Fig.  4(a)  and  (b)  are  in  reasonable  agreement.  Again,  simula¬ 
tions  involving  only  the  polar  polytype  demonstrate  higher  average 
stiffness  and  less  transformation. 

3.4.  Discussion 

Conditions  deemed  sufficient  for  onset  of  amorphization  or 
strength  loss  are  compared  in  Table  4.  For  shock  simulations  (first 
two  rows),  p  is  the  minimum  applied  axial  pressure  at  which  trans¬ 
formation  is  observed.  For  example,  for  simulations  of  the  polar 
polytype  only,  transformation  is  observed  at  p  =  20  GPa  but  not  at 
p  =  1 5  GPa.  Tabulated  results  for  uniform  compression  are  averages 
over  multiple  simulations  discussed  in  Section  3.3.  Prior  analyti¬ 
cal  results  (Clayton,  2012)  involve  elastostatics  and  uniform  strain 
among  grains.  Comparison  of  these  results  with  those  in  the  first 
and  third  rows  indicate  that  stress  dynamics  and  elastic  grain  inter¬ 
actions  promote  instability  relative  to  the  static  and  homogeneous 
case  in  Clayton  (2012).  Experimental  plate  impact  data  are  aver¬ 
ages  over  multiple  shots  (Vogler  et  al.,  2004;  Zhang  et  al.,  2006), 
with  P  corresponding  to  the  average  HEL,  above  which  a  measurable 
strength  loss  is  evident  (Vogler  et  al.,  2004),  though  amorphization 
has  not  been  definitively  proven  to  cause  such  strength  reduc¬ 
tion.  Stress  P  for  ballistic  penetration  experiments  is  estimated  in 
Chen  et  al.  (2003),  with  u  the  projectile  velocity.  For  v  >  900  m/s,  a 
marked  decrease  in  ballistic  penetration  resistance  was  noted,  and 
the  boron  carbide  target  cleaved  into  numerous  fragments  contain¬ 
ing  amorphous  zones. 

Comparing  the  present  results  (first  four  rows)  with  others  in 
Table  4,  satisfactory  agreement  is  found  for  simulations  involving 
the  polar  polytype.  Transformation  is  premature  (i.e.,  p  is  deemed 
too  low)  and  f  too  high  relative  to  experiment,  when  the  CCC  poly¬ 
type  amorphizes  upon  attainment  of  intrinsic  B-instability  criterion 
(2.10).  Thus,  the  present  work  suggests  that  intrinsic  instability  cri¬ 
terion  (2.10)  is  adequate  for  initiation  of  amorphization  in  the  polar 
polytype,  but  not  the  CCC  polytype.  This  result  supports  the  finding 
from  DFT  calculations  (Yan  et  al„  2009)  that  a  buckling  instability 
during  [00  01]  compression  in  the  polar  polytype  leads  to  cross- 
linking  of  the  B  atom  of  the  C-B-C  atom  chain  of  the  rhombohedral 
unit  cell  with  one  of  the  icosahedra  in  the  surrounding  structural 
cage,  followed  by  structure  collapse  and  amorphization.  Since  no  B 
atom  exists  in  the  CCC  poly  type  (i.e.,  its  atomic  backbone  is  C-C-C), 
the  same  mechanism  is  not  possible  in  crystals  of  the  CCC  polytype 
in  B4C. 

The  amorphization  model  implemented  here  is  simple,  with 
no  adjustable  parameters,  and  appears  adequate  for  addressing 
initiation  of  transformation.  However,  some  experiments  (Chen 
et  al.,  2003;  Subhash  et  al.,  2008;  Yan  et  al.,  2009)  report  thin  bands 
(several  nm)  of  glassy  phase,  features  that  cannot  be  fully  resolved 
by  FE  representations  of  polycrystals  with  minimum  element  sizes 
on  the  order  of  tens  to  hundreds  of  nm.  An  alternative  numerical 
technique  is  needed  to  resolve  such  small  features.  Other  nonlinear 
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Fig.  6.  Uniform  compression  of  polar  polytype  at  V/Vo  =  0.95:  (a)  axial  stress,  (b)  shear  stress  and  (c)  order  parameter. 


models  have  addressed  thermal  effects  (Clayton,  2009,  2011b) 
which  could  be  important  in  shear  localization  (Grady,  2011)  but 
are  not  treated  here.  Plastic  slip  and  twinning  have  also  been 
analyzed  for  other  ceramics  in  finite  deformation  models  (Clayton, 
2009,  2011b),  but  (partial)  dislocations  are  thought  to  be  relatively 
immobile  in  typical  B4C  (Clayton,  2012).  Twin  boundaries  are 
scarce  in  conventional  boron  carbide  produced  by  hot  pressing, 
but  are  abundant  in  the  less  common  slip-cast  material  and  affect 
its  stress  fields  (Sano  and  Randow,  2011).  Fracture,  including 
grain  boundary  failure  and  especially  intragranular  (i.e.,  cleavage) 
fracture,  inevitably  occurs  in  ballistic  impact,  and  accompanies 
amorphization  in  diamond  anvil  cell  (DAC)  experiments  (Yan  et  al., 
2009).  Fracture  in  boron  carbide  -  which  could  be  modeled  using 
cohesive  zone  (Clayton,  2005;  Clayton  et  al.,  2012a)  or  continuum 
damage  mechanics  (Clayton,  2006;  Aslan  et  al.,  2011)  techniques 
-  would  presumably  affect  local  stress  states  and  amorphization. 

The  present  limited  results  suggest  that  for  similar  values  of 
applied  or  average  compressive  stress  p,  instability  sometimes 
appears  more  likely  in  shock  simulations  than  uniform  compres¬ 
sion  simulations.  For  example,  transformation  in  the  polar  polytype 
occurs  at  shock  stresses  between  15  and  20GPa,  but  does  not 
occur  in  uniform  compression  until  23  GPa  (second  and  fourth 
rows  of  Table  4).  Statistics  from  more  simulations  involving  more 
microstructures  are  needed  to  verify  this  assertion,  and  the  same 
conclusion  cannot  be  affirmed  from  the  present  results  for  simu¬ 
lations  involving  a  mixture  of  polytypes.  According  to  the  present 
amorphization  model,  for  the  same  microstructure,  differences  in 
transformation  behavior  would  be  due  solely  to  differences  in  local 


stress  state.  In  real  experiments,  entropy  and  temperature  rise 
during  shock  loading  would  tend  to  exceed  that  during  uniform 
isentropic  compression  (Thurston,  1974),  and  this  could  promote 
instability  in  a  way  not  accounted  for  by  the  present  simulations 
that  omit  thermal  effects. 

Several  other  issues  warrant  further  discussion.  Results  demon¬ 
strate  that  the  FE  mesh  resolution  is  adequate  to  resolve 
transformation  zones  and  stress  concentrations,  with  multiple  ele¬ 
ments  tending  to  span  such  regions.  For  example,  contours  in 
Figs.  2,  3,  5  and  6  show  that  stress  concentrations  and  transformed 
regions  are  typically  controlled  by  grain  geometry  rather  than 
local  element  sizes  and  shapes,  though  in  some  cases  regions  of 
amorphized  material  approaching  the  element  size  were  detected. 
Simulations  reported  already  were  conducted  with  default  param¬ 
eters  (Jung,  2010)  for  linear  and  quadratic  artificial  bulk  viscosity 
enabled.  Simulations  were  also  explored  with  viscosity  disabled, 
yielding  results  that  differed  little  from  those  shown  in  Section  3. 
In  shock  problems,  artificial  viscosity  is  typically  used  to  establish  a 
numerical  shock  width  large  enough  to  be  resolved  by  the  grid  size 
(Lew  et  al.,  2002);  in  the  present  simulations,  the  imposed  stress 
wave  is  ramped  and  does  not  localize  to  a  width  equal  to  the  local 
element  size.  Nevertheless,  a  more  comprehensive  study  of  effects 
of  mesh  size/geometry  and  viscosity  may  be  warranted  for  general 
applications  of  the  model.  Extension  of  the  amorphization  model  to 
a  phase  field  approach  (Clayton  and  Knap,  2011a, b),  as  mentioned 
in  Section  2.1,  could  enable  regularization  of  widths  of  transforma¬ 
tion  zones  and  ensure  resolution  of  locally  heterogeneous  features 
of  the  solution  that  might  otherwise  depend  on  mesh  density. 


Table  4 

Sufficient  conditions  for  initiation  of  structure  change  and/or  strength  loss  in  B4C. 


Simulation  or  experiment 

Polytype(s) 

Porp  (GPa) 

v  (m/s) 

V/Vo 

Reference 

Shock  simulation 

CCC  + polar 

15.0 

401 

0.964 

Present  work  Section  3.2 

Shock  simulation 

Polar 

20.0 

607 

0.954 

Present  work  Section  3.2 

Uniform  dynamic  simulation 

CCC  +  polar 

13.1 

- 

0.969 

Present  work  Section  3.3 

Uniform  dynamic  simulation 

Polar 

23.0 

- 

0.948 

Present  work  Section  3.3 

Nonlinear  elasticity  analysis 

CCC  +  polar 

18.1 

- 

0.960 

Clayton  (2012) 

DFT  single  crystal  calculation 

Polar 

18.9 

- 

0.930 

Yan  et  al.  (2009) 

Plate  impact  experiment 

Not  reported 

16.5 

600 

0.944 

Vogler  et  al.  (2004) 

Plate  impact  experiment 

Not  reported 

19.5 

530 

0.965 

Zhang  et  al.  (2006) 

Ballistic  penetration  experiment 

Not  reported 

23.3 

900 

- 

Chen  et  al.  (2003) 
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4.  Conclusions 

A  geometrically  nonlinear  theory  for  anisotropic  elasticity  and 
stress-induced  amorphization  has  been  implemented  for  boron 
carbide  and  exercised  in  dynamic  compression  simulations  of 
polycrystals.  Results  for  stress-strain  behavior,  initiation  of  trans¬ 
formation  in  the  polar  polytype  associated  with  loss  of  intrinsic 
stability,  and  strength/stiffness  loss  compare  favorably  with  atomic 
theory  and  experimental  observations. 


that  definition  (A.7)  is  exact  when  S2  E->0  (Clayton,  2012).  This 
criterion  agrees  exactly  with  that  for  classical  stability  under 
hydrostatic  loading  (Milstein  and  Hill,  1979).  Correlations  are 
likely,  but  not  always  necessary,  among  failure  of  intrinsic  stability 
criteria  in  perfect  crystals  and  the  onset  of  various  structural 
changes,  e.g.,  melting,  cavitation,  slip,  twinning,  or  solid-solid 
phase  transitions  (Born,  1940;  Hill,  1975;  Parry,  1980b;  Wang 
et  al„  1995;  Morris  and  Krenn,  2000). 

References 


Appendix  A.  Intrinsic  stability 


Consider  an  elastic  material  with  strain  energy  per  unit  refer¬ 
ence  volume  <F(  q )  undergoing  homogeneous  deformation.  Let  S'!* 
denote  an  increment  in  strain  energy,  and  Sw  a  first-order  incre¬ 
ment  in  external  work  at  fixed  p,  where  q  and  p  are  vector 
representations  of  generalized  coordinates  (i.e.,  deformation)  and 
conjugate  force: 


d&  1  3 V 

Sq,=  °.Sq+-Sq.-^- 
dq  2  dqdq 


Sq  +  ■ 


Sw  =  p-Sq. 


(A.l ) 


Classical  stability  requires  >  Sw  (Hill,  1975;  Parry,  1978).  At 
fixed  load  (i.e.,  stationary  virtual  stress  as  in  (A.l ))  with  first-order 
equilibrium  constraint  p  =  9<F/9  q  and  to  second  order  in  S  q, 
“intrinsic”  stability  criteria  emerge: 


d2<P 

8q^-8q>0 

dqdq 


(V<7  =f=  0) 


det 


92<F 

dqdq 


>  0. 


(A.2) 


Such  criteria  (which  amount  to  local  convexity  conditions),  while 
dependent  only  on  the  local  state,  depend  non-uniquely  on  choice 
q  and  correlate  with  classical  stability  only  in  certain  environments, 
e.g.,  load  control  via  p.  Let  8  e  =  F~T8  E  F^1 ,  where  8  E  is  an  incre¬ 
ment  in  Green  strain.  Incremental  external  work  per  unit  reference 
volume  is 


8w  =  S-.8E  =  P:8F=J<r:8e.  (A.3) 

Born’s  mechanical  stability  criterion  (Born,  1940)  involves  the 
choices  and  q-t  £  and  p-»  S  =  9*F/9f,  with  spatial  (c)  and  refer¬ 
ential  (C)  tangent  elastic  moduli  related  by  Cy ki  =J_1  FiiFjjFki<FiLQji<L- 

32<F 

8S:SE  =  SE:  T-r-  :SE  =  8E:C:SE 
dEdE 

=  J8e  :  c  :  Se  >  0  <s>  det[C]  >  0  det[c]  >  0.  (A.4) 


For  intrinsic  stability  consistent  with  classical  stability  under  all- 
around  dead  loading  (Hill,  1975;  Ogden,  1984),  F  and  p-> 
P  =  3<F/3F: 

32<P 

8P:8F  =  8F:  :  8F  =  SF  :  A  :  8F  >  0  <s>  det[A]  >  0; 

oFoF 

A-iJkL  =  FuFkKCyKL  +  SjLS  jk-  (A.5) 

For  intrinsic  stability  consistent  with  classical  stability  under  con¬ 
trolled  Cauchy  stress  (Wang  et  al„  1995;  Morris  and  Krenn,  2000), 
p-s*  cr  and  q-s*  e: 


Sa  :  Se  = 


:  Se  =  Se  :  B  :  Se  >  0  det[B]  >  0; 


(A.6) 


Bijkl  =  ci]kl  +  2  (aikfyl  +  ail^jk  +  VjAk  +  ajAl  ~  &ij8kl  ~  (A.7) 


Incremental  spatial  modulus  B  was  derived  45  years  ago  by 
Wallace  (1967)  but  has  only  been  used  more  recently  in  general 
stability  criteria.  The  first  two  equalities  in  (A.6)  strictly  apply  only 
when  F  is  symmetric,  but  the  inequalities  apply  regardless  since 
B  is  rotationally  invariant.  In  the  former  case,  it  can  be  proven 
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